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Abstract 

The grid cells (GCs) of the medial entorhinal cortex (MEG) and place cells (PGs) of the hippocampus are 
assumed to be the key elements of the brain network for the metric representation of space. Existing the¬ 
oretical models of GC network rely on specific hypotheses about the network connectivity patterns. How 
these patterns could be formed during the network development is not fully understood. It was previously 
suggested, within the feedforward network models, that activity of PCs could provide the basis for devel¬ 
opment of GC-like activity patterns. Supporting this hypothesis is the finding that PC activity remains 
spatially stable after disruption of the GC firing patterns and that the grid fields almost disappear when 
hippocampal cells are deactivated. Here a new theoretical model of this type is proposed, allowing for grid 
fields formation due to synaptic plasticity in synapses connecting PCs in hippocampus with neurons in MEG. 
Learning of the hexagonally symmetric fields in this model occurs due to complex action of several simple 
biologicaly motivated synaptic plasticity rules. These rules include associative synaptic plasticity rules simi¬ 
lar to BCM rule, and homeostatic plasticity rules that constrain synaptic weigths. In contrast to previously 
described models, a short-term navigational experience in a novel environment is sufficient for the network 
to learn GC activity patterns. We suggest that learning on the basis of simple and biologically plausible 
associative synaptic plasticity rules could contribute to the formation of grid helds in early development and 
to maintenance of normal GCs activity patterns in the familiar environments. 


1. Introduction 


The grid cells (GCs) in the medial entorhinal cortex (MEG) of mammals are the neurons whose firing 
activity during animal navigation is concentrated near the centers of hexagonally symmetric grid covering 
the environment (Hafting et al. 2005 Eyhn et al., 2007| Stensola et al. 2012 Moser et al. 2008, 2015). The 
network of GCs is considered to be a component of the brain system for metric representation of space and 
navigation (Moser et al., 2008 Hartley et al. 2014[ [Moser et a~ 20151. Patterns of activity, connectivity 
and development of GC networks have been extensively investigated, and a number of theoretical models 


of GC network was developed during last decade (see ( 

Zilli 

20121 

and (( 

Hoser et al. 

2014 

1 for review and 

(Pilly & Grossberg 

20131 Castro & Aguiar 

2014 IWidloski & Fiete 

2014 

Bush & BurgessI 2014| Hasselmo 

& Shay 

2014 

1 for recent models). These models differ significantly in the way they represent current 


position of the animal. For example, models based on continuous attractor networks (CAN) (Conklin & 
Eliasmith 2005 McNaughton et ^ 2006 O’Keefe & Burgess, 2005 Fuhs fc Touretzk^ |2006|) , interf erence 


of slow oscillations (Blair et al. 2007 Burgess et al. 2007 Gaussier et al. 2007 Giocomo et al. 


and feedforward networks (Kropff & Treves 2008 Castro & Aguiar 2014 Pilly & Grossberg 20131 exist. 


20071 


Depending on the nature of input signals these models could be divided into path-integrating models and 
models that utilize location-specific combinations of input signals. Currently, none of the existing GC 
network models could be considered as sufficiently complete and supported by the experimental data. 

Most of the GC network models could be assigned to the path-integrating type. They imply the emer¬ 
gence of a network, controlled by the signals about animal’s speed and direction of movement, in order to 


provide integration of movements along the path (Zilli, 2012 Moser et al. 20141. Models of this type have a 


number of restrictions. For example, path integrating models use only a small part of sensory input signals 
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from those which are received by EC and could be the source of navigational information. Other set of 
problems arise because path integration process in these models is assumed to be continuous, and gradually 


accumulating integration errors should constantly be adjusted by some additional mechanism (Hardcastle 


et al. 

2015 

Cheung et al. 

2012 


Most of the existing models of path-integrating type rely on relatively complex predefined patterns of 
network connectivity and do not address the problem of these patterns development. For example, many 


path-integrating models, which use CAN for animal position representation (Conklin & Eliasmith 2005 


McNaughton et al. 2006 O’Keefe & Burgess, 2005 Fuhs & Touretzky 2006), are based on hypotheses that 


1) two-dimentional neuronal network with toroidal topology of synaptic connections exists in the entorhinal 
cortex, 2) strength and symmetry of these connections are high enough for a stable localized bump or 
a periodic set of bumps of activity to exist in the network, and 3) additional groups of neurons exist 
with activity, modulated by speed and direction of animal movement, and with asymmetric connectivity 
within attractor network layer. As a result, the bump of activity in 2D layer of neurons could move in a 
coordinated manner with animal movement. Several possible scenarios of development of path integrating 
CAN were proposed (Stringer et al., 2002 Hahnloser, 2003 McNaughton et al. 2006 Widloski & Fiete, 2014). 
The problem of development of most of connectivity patterns, assumed by hypotheses l)-3), was directly 
addressed in a recent model of Widloski and Fiete (Widloski & Fiete 2014). In their work, formation of 
path-integrating 2D network of GCs under control of both location-specific and velocity-modulated input 
signals was demonstrated. This network was obtained as a result of self-organization process driven by 
spike-time dependent plasticity (STDP) in the lateral synaptic connections in the network of excitatory and 
inhibitory neurons. 

Another type of the GC models suggests that the feedforward network connecting neurons with spatially 
modulated patterns of activity, such as hippocampal place cells (PCs) or MEC/parasubicular stripe cells. 

Si et al. | [M^[PiII 71 F 


to GCs is formed due to associative synaptic plasticity. (Kropff & Treves, 2008 


Grossberg 2013 Castro & Aguiar 2014). In several models of this type, formation of grid fields in a given 


environment is associated with formation of synaptic connections to GCs selectively from those PCs whose 


firing fields are located in a nodes of hexagonal lattice (Kropff & Treves 2008 Castro & Aguiar 2014). 


This connectivity pattern is supported by the finding that PC activity remains spatially stable after the 


disruption of GC firing patterns (Koenig et al. 2011 Hales et al. 2014 Bush et al. 2014) and that grid 
patterns almost disappear when hippocampal cells are deactivated (Bonnevie et al., 2013). Development 


of the early PCs before GCs formation also supports the role of PCs as spatial information providers to 


GCs.(Langston et al. 2010 Wills et al. 2010) 


The limitation of models based on PC-GC feedforward network self-organization (Kropff & Treves 2008 


Castro & Aguiar, 2014) is that they do not explain explicitly how normal pattern of GC activity could 
persist during animal’s navigation in darkness (Hafting et al. 2005). Other GC network property which 


seems to be difficult to explain within PC-GC feedforward network hypothesis is the fact that in a novel 
environment relative phases of grid fields in the pairs of GCs do not change, while place fields rebuilt 


completely and their relative phases in PC pairs before and after remapping are not correlated (Fyhn et al. 


2007). As it was described in the work of Kropff and Treves (Kropff & Treves 2008), alignment of grid 


fields orientations in the network and correlated remapping of grid fields in different environments could be, 
in principle, achieved if GCs activity before grid fields formation is modulated by head direction cells and if 
the strength of horizontal excitatory synaptic connections, which form between pairs of GCs in 2D layer of 
GCs, is modulated by the degree of similarity of their prefered head directions and by their relative distance 
in 2D layer. In the work of Si, Kroppf and Treves (Si et al. 2012) it was shown that coordinated remapping 
of grid helds between two environments indeed can be observed in the model if the sets of presynaptic PCs 
active in these environments are weakly overlaping. Modular organization of GC networks, and the fact that 
formation of new grid fields in a novel environment usually do not require prolonged sensory experience in 
this environment, are the other issues that require explanation within the framework of feedforward network 
models (Barry et al. 2007 Stensola et al.) 2012 Barry et al. 2012). 

Models based on PC-GC feedforward network self-organization rely on specific plasticity rules and differ 
by their characteristics. For example, the model of Kropff and Treves (Kropff & Treves, 2008 Si et al. 


2012) assumes that the mechanism of GCs firing rate adaptation besides Hebbian plasticity rules in PC-GC 
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synapses is important for grid fields development. The model of Castro and Aguiar (Castro & Aguiar 20141 


assumes the existence of a specific complex rule linking synaptic states and animal’s position. At the same 
time, it is well-known that the patterns with hexagonal symmetry are one of the most ubiquitous types of 


patterns that emerge as a result of symmetry-breaking bifurcations in many physical systems (Golubitsky 


et al. 19881, in particular in the models of self-organization based on Hebbian plasticity in neural networks 
described by the neural fields equations (Ermentrout & Cowan, 1979 Bressloff et ah, 20011. Thus, it is 
interesting to explore other possible mechanisms of grid fields development, with a more general types of 
plasticity rules involved. 

The goal of this work is to describe a class of models in which the grid fields are formed as a result of 
plasticity of PC-GC synapses or synapses from sensory neurons to GC, and to demonstrate a number of 
possible grid fields formation scenarios in these models. 

The main difference of the model proposed in this paper from the previously proposed models is that it uses 
only simplest rules of synaptic plasticity, similar to those supported by experimental data, and does not 
require assumptions about the firing rate adaptation. 

The most important assumption made in this work is that dependence of associative plasticity rate of 
PC-GC synapses on the frequency of presynaptic action potentials has a minimum within a typical range 
of PCs activity during animal’s exploration of environment. This assumption is supported by a number 
of experimental works demonstrating that transition from long term depression (LTD) at low pre- and 


postsynaptic firing rates to long term potentiation (LTP) at high rates could occur in hippocampal ( 

Bear & 

Malenka 

1994 

Artola & Singer 

1993 

Dudek & Bear, 1992 

Mayford et al. 

1995 

Wang & Wagner| 

1999 

1 and 

entorhinal cortex synapses ( 

Solger et al. 

2004 

Zhou et al. 

2005 

Deng & Lei 

20G 

71. The second assumption 


is that some components of associative or homeostatic synaptic plasticity depend nonlinearly on synaptic 
weights, and that the curvature of this dependence is positive. This assumption seem to be physiologically 
relevant since it simply results, for example, from lower-bounding of exitatory synapses weights. 

Learning of grid fields in the proposed model could be very fast and does not require a lot of exploratory 
experience in a novel environment. In addition, the proposed model is compatible with CAN models of GC 
network and could serve as a basis for the path integration error correction, and for the self-organization of 
activity bump position control in them. 


2. Theory 

2.0.1. Description of the model 

In this work we develop a hypothesis that grid cells recognize sets of sensory signals, associated with 
location of an animal in certain places of environment, with the help of clusters of synapses from hippocampal 
place cells to grid cell or from neurons of different sensory systems to grid cell. Clusters are composed of 
synapses with correlated presynaptic activity and form as a result of Hebbian associative plasticity. In 
addition, under the action of both associative synaptic plasticity and homeostatic plasticity, which limits 
weights of the individual synapses or groups of synapses, separation of clusters occurs in such a way that 
different clusters encode maximally different sets of local sensor signals. 

As a possible mathematical model of such a system we consider a network of neurons with linear activation 
rule 

j/t=wfxt (1) 

where xj is a vector of firing rates of presynaptic neurons at time moment t, w* vector of PC-GC synaptic 
weights and yt - firing rate of postsynaptic neuron. 

We will consider several types of synaptic plasticity rules which will lead to different scenarios of grid 
fields self-organization. Let us start from the general form of equation describing synaptic weights change 
as a function of rates of presynaptic and postsynaptic action potentials: 


Wit = V {xit, yt) + C {yt, wt) Wit + P{w^t,Xit) 


( 2 ) 
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where ( and P terms describe global homeostatic plasticity (restricting sum of weigths of all synapses) 
and local homeostatic plasticity (restricting weights of individual synapses independently of each other). 
The t] term describes associative plasticity corresponding to LTP and LTD. For simplicity, we approximate 
homeostatic plasticity terms by polynomials: 


Ciyt) = -cryt or ^(wt) = -cr^ 


wit, 


P{wit,Xit) = A{xit)wit + B{xit)wlf., 


( 3 ) 

( 4 ) 


where A(0) = 0, 5(0) = 0 and a is some positive number. Here, we assume that local homeostatic plasticity 
rate is relatively low, so that third and higher order terms of its Taylor series expansion could be neglected. 

From experiments it is known that transition from LTD at low pre- and postsynaptic firing rates to LTP 
at high rates could occur in hippocampal (Bear & Malenka, 1994 Artola & Singer 1993[ Dudek & Bear 


1992 Mayford et al., 1995 Wang & Wagner 19991 and entorhinal cortex synapses (Solger et al. | |2004[|Zhou| 


et al. 2005 Deng & Lei 20071. To describe this phenomenon, polynomial approximation of r; (ccit, y*) should 


contain terms up to at least third order: 

V {Xit, yt) = -Tl-Xityt + m+Xityl + V+xftyt 


( 5 ) 


where coefficients ? 7 _ (rate of LTD), 77 + (rate of LTP linear component) and 772 + (nonliner LTP rate) are 
positive. Without last term this equation is similar to the Bienenstock-Cooper-Munro (BCM ) rule 


iienen- 


stock et al., 1982 Cooper fc Bear[ 2012), phenomenological synaptic plasticity rule predicted theoreticaly 
to account for selectivity of primary sensory cortex neurons. Similar plasticity rules were obtained by av¬ 
eraging of plastic changes that occur in response to a Poisson stream of pre- and postsynaptic spikes in 


detailed phenomenological models of STDP (Gjorgjieva et al. 2011 Mayr & Partzsch 2010 Bush et al. 


2010). Other possible interpretations of Eq.([5) is that its first and second terms are equivalent to linear and 


quadratic terms of the Taylor series expansion of anti-Hebbian plasticity rate dependence on postsynaptic 
activity. This could occur, for example, if neuron has nonlinear activation function S{z) 


yt = S{zt) = S{wjxt) 
rj [xit, zt) = -TloXityt = -y-XitZt + r] 2 +Xitzl + o{zf) 


( 6 ) 

( 7 ) 


Here yt is the output firing rate and Zt is the average change in postsynaptic membrane potential evoked by 
synaptic activity. 

Last term in Eq.([^ could be a result of nonlinearity of LTD rate dependence on presynaptic activity. As we 
will see later, relatively small nonlinearity of this dependence could be sufficient for grid fields formation. 
There are experimental results which support assumption that this term indeed could be required for precise 
description of STDP at least at low rates of postsynaptic activity. Eor example, it was shown that spike 


time dependent LTP could be observed in absence of postsynaptic spiking (Bi & Poo, 1998 Nevian & 


Sakmann, 2006). Probably in this case, at low postsynaptic activity levels, the rate of synaptic changes 


could approximately linearly depend on yt and nonlinearly on Xit- Other motivation for introducing terms 
with differences in dependence of LTP component of plasticity on presynaptic activity with respect to LTD 
is experimental results suggesting that LTP and LTD could be selectively supported by different sets of 


postsynaptic receptors and voltage-gated channels (Nevian & Sakmann 2006 Solger et al., 2004). But, as 


we will see later, LTP component is not required for development of grid fields. 


2.0.2. Analysis of weakly nonlinear version of the single neuron model 

To analyze potential mechanism of grid fields formation in single neuron model with Hebbian and anti- 
Hebbian synaptic plasticity described by Eqs. (H-® it is convenient to consider a variant of this model 
which is weakly nonlinear with respect to synaptic weights. Namely, we will assume that in Eqs.(|^-([^ all 
the terms that are nonlinear with respect to are small, except for those which are equal for all synapses 
( cr{yt,wt)wt ). 
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If learning rate is slow with respect to the exploration speed, learning equations Eqs. ©-(Hi could be inte¬ 
grated over t : 


Wi = -Ty-Rij-w -f 77+R2iW -h ? 72 +w'^Qj’w -1- a (w) Wi + P{wi) 

Here Ri,R 2 and Qj are time averages of functions of input signals: 

(8) 

R-i = 

t 

(9) 

t 

(10) 

Q* = 

t 

(11) 

If the animal explores environment for a long time, and visits all places with equal probability, or if learning 
rates are normalized according to the extent of place familiarity, then the matrix Ri could be expressed as 
a function of presynaptic PCs firing fields Fi{r) : 

Ruj ~ J Fi{r)Fj{r)dr° 

(12) 


We assume that all place fields are circulary symmetric within some distance from their centers and are 
approximately equal: 


Pi{^) « /(Ik-ri 


(13) 


For approximately Gaussian place fields -Fi(r) 
analytical expression: 


Fexp(—^ |jr — , cross-correlation Ri has simple 


R 




1 
<7 

V, - r 


Analogously for R 2 and Q : 


R 


'2ij 


j Fi{rfFj{r)dr^ = 2 tt^F^ exp ll^i - 

Qijk^ [ Fi{r)Fj{r)Fk{r)dr^ = {RiijRukRijky 


(14) 

(15) 

(16) 


It is known that development of periodic patterns with different symmetries may be observed in systems 
described by equations similar to Eq.([^ as a result of symmetry-breaking bifurcations. The symmetry 
of periodic pattern, which begins to grow in this and more general systems after the bifurcation, can be 


predicted using neural fields approximation and methods of weakly nonlinear analysis (Bressloff 2005 


Ermentrout 1991 Tass 1997). The stability of stationary periodic patterns with respect to their periodic 


perturbations could be investigated using methods of symmetry-breaking bifurcation analysis for neural 
fields equations (Golubitsky et ah, 1988 Bressloff et al. 2001). 

Here we only briefly discuss, using standard approach, the mechanism of formation of spatialy periodic 
sensory fields with hexagonal symmetry in the system Eq.([^. First, note that Lyapunov function exists for 
Eq.(|8|: 


E = rj- -w^Ri-w — -r74.’W^R2W — \r] 2 + f {w'^Xrf dr^ + ja {w^wf f P{w,)dw, 


(17) 


As a consequence, all solutions of Eq.( [^ s hould converge to a steady-state solution that provides a local 
minimum of the Lyapunov function Eq.(17). 
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This local minimum always exists since the length of weight vector L^, = Vw’^w will decrease for all 
weigth vectors with the length larger then some limit. For example, for Eq.Q without third term (i.e. in 
the case 772 + = 0 ) this follows from the set of inequalities: 


dt 



— w^w = w^RoW — tr (w^w) + w^P(w ) ^ ^^'^rnax 

^ ^maxL^ju ~ + BNL^ 


(18) 


where Xmax is the maximal eigenvalue of the matrix Rq = —? 7 _Ri +? 7 +R 2 +eiE, N- the numbet of elements 
in the vector w and Wmax - the element of w with a maximal absolute value. It can be seen that < 0 
for all Lyj larger than the largest root of the polinomial in the right side of Eq.(18). 

Futher analysis will show that the solution of Eq. Q will always approach grid-like pattern at least for 
some sets of parameters. We assume for simplicity that the domain for place fields spatial locations is a 
twisted torus. In other words, each place field F(r) has periodic structure 


F(r) = F(r + mviL-I-n 2 V 2 L) (19) 

where ni,n 2 G Z, Vi = (1,0) and V 2 = (l,^) basis vectors of the hexagonal lattice of size L. In 
this case, first three terms in equation Eq.® are convolutions. As follows from the convolution theorem, 
the equation for the spatial Fourier decomposition of weights distribution w as a function of place fields 
positions r has simple form: 


= E' 


Wr = > , Q!ke 
keA 


..ikr 


ttk = 


—77_Rl|k|Q;k + ??+R 2 |k|Q:k + rj2+ E] Ql|k|,|k'|,|k-k'|Q:k'Q:k-k' 

k' 

-cr ^ Ickk' l^ttk -t- Pi (a, k) 

k' 

Pi {a, k) = EiCkk -I- £2 Clk'lAk-k' 


( 20 ) 

( 21 ) 

( 22 ) 


k' 


Ikl 


ll '1^ cr„ 

|k| 


k-|-k'|^^ for place fields described by Eq.(14), ei « / A{F{r)) dr^, €2 


where wavevector k is defined on the hexagonal lattice A which is dual to the lattice A, Qi|k| |k'| |k-k'| = 
FqCxp 

jB{F{r))dr^. 

Let us initially consider evolution of wave amplitude, ctk , starting from small positive initial values 
of synaptic weights. If ei is positive, then there is some k for which —77_Ri|k|0!k + ? 7 -i-R 2 |k|Q:k + eicck is 


positive as well. Since for small ak influence of all high order terms in Eq.(21) can be neglected, then the 
most rapidly growing wave amplitudes are those located on the circle |k| = kc, where kc is the amplitude 
of wavevector for which expression —77_R]^|k| -I- 77+R2|k| is maximal. As a result of exponential growth of 
wave amplitudes, to the end of this ’linear’ phase only wave amplitudes with |k| = kc will be significantly 
different from zero, but amplitudes of all these waves grow with approximately equal speed. In other words, 
the ratios of amplitudes of waves with |k| = kc do not change with time if 772 + = 0 and £2 = 0. 

The influence of second order terms gradually increase with growth of ak, but if 772 + and £2 are small, 

then most of the wave power will still be concentrated within the narrow ring Ro,|k| ~ R-o,fec ^ ^2 > where 
Ro.|k| = ~^-Rl|k| + ?7-l-R2|k| + £l- 

The action of second order terms on the wave growth generally results in gradual increase of differences 
in speed of growth of different wave groups, so that only small-sized group of waves will dominate at the 
end. Indeed, let us consider a group of waves with wavevectors lying on the circle |k| = kc- For the first- 
order approximation, the contribution of all |k| ^ kc to the second order (with respect to ak ) terms of the 


equation Eq.(21) can be neglected. 
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For each wavevector k lying on the circle |k| = kc the re is only one triple of wavevectors k, k^, k-k^ with 
|k| = Ik'l = |k-k'| = kc- Consequently, the equation Eq.(21) could be simplified: 


dk — |^R-0,|k| ~ ^ lo^k'I J CUc + 2Q2_|k|,|k+|,|k-|C^k+Q:k- 

Q2,|k|.|k'|,|k-k'| = ^2 + ^2+Qi,|k|.|k'|,|k-k'| 


(23) 

(24) 


where k+ and k— are wavevectors obtained by rotation of wavevector k by an angle ^ and — ^ , respectively. 
Since synaptic weights are real-valued, then amplitude of wave with wavevector k is equal to the complex 
conjugate of the amplitude of wave with opposite wavevector: a_k = dk- We can observe that differences 
between absolute values of amplitudes of waves, having wavevectors rotated by the angle ^ with respect to 
each other, will decrease with time. Indeed, let = |a_k|e®'^'' and let the vectors k2-|- and k2— correspond 


to wavevector k rotated by angles ^ and , respectively. Then from Eq.(23) it follows that: 


d 

dt 


l«k| = (^Ro.lki - l«kf ^ |ak| + 2Q2,|kMk2+|.|k2-| I«k 2 +| |ak 2 -| i?e 

<^k = -2Q2,|k|.|k2+Mk2-| l«kr' |ak2+| |ak2-|/m{e*(^>‘+^-++^--)} 


As a result, at t —> oo the sum of phases ((<k + </'k 2 + + 


'k2- 


0 and Eqs.([^ converge to: 


dt 


l^kl I Ro,|k| ^ E ] Ic^kl + 2Q 2,|k|,|k2-|-|,|k2-| |Q:k2+l I®k2-| j 


(25) 

(26) 

(27) 


k' 


This equation describes growth of absolute values of amplitudes, that will continue at least while Ro,|k| ~ 
a |ak' 1^ > 0. From some time moment the last inequality will be broken for all t > Td and differences 
between absolute values of wave amplitudes for triplets of wavevectors (k2, k2—, k2-|-) will decrease. Indeed, 


from Eq.(27l we have: 


|(K|-|ak2-H|): 



• a 


E 

k' 


■ 2Q2,|k|,|k2+|,|k2-| I«k2-| 1 (lokl — |ak2+|) , (28) 


Since expression in left brackets on the right side of Eq.(28) is always negative for t > Td then |ak| ~ 
|Q:k 2 -i-l —0 with time. As a result, in a steady-state all nonzero wave amplitudes should be equal. 

In general case, the number of wavevectors defined on the hexagonal lattice A and lying on the circle 
||k|| = kc, is an arbitrary natural number divisible by six. But only the configuration with six nonzero waves 
whose wavevectors compose equilateral hexagon is stable. 

To show this, let us consider evolution of ratios of wave amplitudes Uk = where aki is the amplitude 
of the wave having maximal amplitude at t —> oo : 


Ak — Q!kl ( 2CI2, kc,kc ,kc + ‘^^ 2 ,kc,kc,kcki(.+ kk-) 


As follows from Eqs.(25 -^81, at t —)■ oo this expression approaches: 


Ak — ‘^^ 2 ,kc,kc,kc^^^kit- {kit 1) 


(29) 


(30) 


This equation shows that if at some moment of time |/Xk| becomes less than I, for example due to some small 
perturbation, then if k do not belong to the group of wavevectors obtained by rotations of kl by angles ^n, 
n G Z then |/rk| will be negative until has reached zero. 


We can conclude that for sufficiently small 772 + and 62 the solution of Eq.(21) , which starts from small 


positive weights, will approach a pattern that is close to the hexagonally symmetric sum of three plane waves: 


Wr = 


aRe I -f gi(r-ro)k 2 +^ | 


(31) 








The same proposition is true for any other starting distribution of synaptic weights. Using first-order 
perturbation approximation with respect to 772 + and £2 it could be shown that steady-state solution of 
Eq.(21) can be approximated by steady-state solution of Eq.(23) and that hexagonally symmetric sum of 
plane waves gives the only stable solution of Eq.(23l and Eq.(211 for sufficiently small 772 + and £2 

If the size of environment is large with respect to the size of place fields, and coefficients near second 
order terms are relatively large, then deviations from the twisted torus assumption will not significantly 
disturb general form of solution of equation Eq.(211, since in this case triples of wavectors lying sufficiently 
k-ku = kc always could be found for each wavector k lying within the ring of large 


k'l = 


As a result, Eqs.(23 - 31) will be approximately correct if Q 2 |i(.| |k'| |k-k'| i® sufficiently 


close to |k| = 
amplitude waves, 
large. 

The size of grid fields formed in the system described by Eq. Q is determined by the critical wavelength 
Ic = 2 tt/ kc- It has a simple physiological meaning: it is equal to a twice distance from the center of place field 
to the position of the animal where rate of associative synaptic plasticity is minimal. For example, in the 
case when 772 -I- = 0, 4 is defined by the distance to a place where the presynaptic firing rate is equal to the 
position of the minimum of associative synaptic plasticity rate dependence on presynaptic firing frequency: 
Ic = 2a7’(777i77r||r|| (—77_i?i(||r|j) -|- r]+R 2 {\\r\\)). In this case grid size depends only on the size and shape of 
place fields and on 77 -/ 77 + ratio. 


2.0.3. Alignment of grid field orientations in models of GCs network 

A network of GCs interconnected by a stable synaptic connections with 2D topology and described by 
Eqs. w could be modeled using activation equations: 


Vgt = -T/gt -f JgYt -bwJjXt 

T ^ 


(32) 


where Jggi = J(||sg — Sgi||) is connectivity matrix of the network and Sg is a position of neuron g in the 
2D layer of GCs, r is membrane time constant of GC, which is small with respect to the time constants of 
learning. If activity of PCs changes slowly with respect to r, and the strength of excitatory synapses is not 
too high, then ijgt « 0 and: 


wxt (33) 

where w is the connectivity matrix for PC-GC synapses. This approximation assumes that all eigenvalues 
of the matrix C = (+E — J) are negative. The learning equations for GC network could be obtained 
from Eqs.Q-0. After setting C(wt) = X)g,fc w^gfct and P{wgit,Xit) = BH(xit - 0)wA^ in Eqs.0-Q we 
obtain: 


Wg,t = -r]-xnygt + m+XgitVgt + O+x^tygi - '^I'ktWgzt + Bw%^H {xgn - O) 

g',k 


(34) 


where H is the Heaviside function. Integration of this equation over time gives, similarly to Eq.Q: 

Wi = — 77 _CwRii -I- 77 +CwR 2 i -I- 772 +Tr |CwQjW^C’^| -|- a (w) Wi -|- P{wi) (35) 

For the simplicity of futher analysis let us consider 772 + = 0 case. Using the same assumptions as for 
derivation of Eqs.(21) after Fourier decomposition of Eq.(35) we obtain: 


keA 


Pglu — C^0|k|/^k ^2 ^ ^ /^^k^/^^k-k' ^ ^ ^ |/^5'k' | Pgk 


g'M' 


(36) 

(37) 
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If we assume that similarly to the layer of PCs, the GCs positions in 2D layer compose a hexagonal lattice 
As, and the domain for the Eq.(35) in the GCs layer is a twisted torus, then for the Fourier decomposition 
of function j3g\^ over the dual lattice As we have: 


Q \ ^ ^ 

Psk — / ^ 

m^As 


^mk — R,o|k|Q^mk “t“ 62 


E' 

k' 


'k'^m-m'k-k' ^ 


E 

m' .k' 


^mk 


(38) 

(39) 


The simplest case for the analysis of the pattern formation in system Eq([M| is the case when connectivity 
matrix J has a maximum at the main diagonal, or in other words function J(||sg — Sgi||) has maximum at 
II Sg — Sgill = 0 . Note that this condition imposes no restriction on the sign of synaptic weights and thus 
such connectivity pattern is possible for both purely excitatory and inhibitory GC-GC synaptic connections. 
In the case of this connectivity pattern, the function C|m| will reach a maximum at m = 0 , and as a result, 
the most rapidly growing wave amplitudes are aok for the wavevector k located on the circle |k| = kc- Using 


arguments similar to those used for derivation of Eqs.(23 - 31) we could obtain a set of similar equations: 


bok 


CnR, 


OtVQ,|k| 


-E 


aok + 2e2bok2+bok2- 


with a solution, approaching: 


= aRe 


{(' 


; _|_ gi(r-ro)k2- gi(r 


-ro)k2+^ I 


(40) 


(41) 


This means that if the strength of GC-GG connections is sufficient, then after learning grid fields of all GCs 
will have the same orientations and phases. It could be hypothesized that in the case of less symmetric or 
more noisy connectivity matrix J the drift of phases could be observed at scales larger than a characteristic 
length of connection in the GCs layer. 


3. Results 

To confirm previous analysis we have conducted a series of computational experiments. First, we have 
found that formation of hexagonally symmetric fields can be observed for a broad ranges of the second order 
terms coefficients 772 + and £2 in Eq.Q. Typical results of computational experiments are shown in Fig. 

Linear neuron with 961 synapses was modeled. Each synapse received spatially modulated input from a 
PC with Gaussian-shaped firing field centered on a vertex of 31x31 hexagonal lattice embedded in diamond¬ 
shaped environment. It can be seen that during formation of grid fields (Fig. upper row) wave amplitudes 
(Fig.j^A, middle row) rapidly concentrate on a narrow ring, which corresponds to the critical wavelength, 
and then selection of single group of six waves of equal amplitudes and with wavevectors composing a 
hexagon occurs. Corresponding autocorrelation functions of grid fields are shown in Fig. (lower row). 
Dependence of gridness score on simulation time step is shown in Fig. Pearson’s correlation coefficient 
of autocorrelation function with the same function after rotation was calculated. It is shown as a function 
of rotation angle in Fig. Gridness score was estimated as a difference between the mean value of this 
correlation coefficient calculated for rotation angles 60° and 120° , and its mean value for angles 30°,90° 


and 150° ( 

Sargolini et al.l|2006 

Si et al. 

2012 

1. 

Model parameters (see Eq. ( 

21)) used for simulations shown in Fig.[^were the following: cr = 1, £2 = 50, 


r] 2 + = 0, El — 0, 77 + = 1, ? 7 _ = I and Ur = 0.08L. With these parameters maximal contribution of second 
order (with respect to synaptic weights) terms to the weight change rate was less than 4% of first order 
terms contribution. Formation of grid fields with high gridness score (more than I) could be observed for a 
broad range of values of parameter £2 - from 0.001 to 100. At larger £2 values unimodal fields usually form 
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Figure 1: Self-organization of grid fields in a simplified model of neuron with associative plasticity. Example of development 
of hexagonally symmetric fields in the model described by Eq. § in the twisted torus domain. A. Evolution of grid fields 
(upper row, spatial distribution of firing activity in central part of environment is shown, height of the square is equal to the 
height of diamond-shaped environment) wave amplitudes (i.e. distribution of the absolute values of variables from Eq. l |20[ l 
in wavevectors space k , middle row) and autocorrelation functions of grid fields (lower row). B. Dependence of gridness score 
on simulation iteration number. C. Pearson’s correlation coefficient of autocorrelation function with its rotation transform is 
shown as a function of rotation angle (([Sargolini et al)||2006| |Si et ah]|2012[l). 


instead of grid fields. As a result, those contribution of second order terms to the weight change rate for 
which grid field patterns formation can be observed generally do not exceed 10%. Similar results can be 
observed in the case when only BCM-like second order term is nonzero {r] 2 + > 0, £2 = 0). 

The formation of grid fields by the mechanism described above could be observed in the model Eq. ([^ for 
the cases when restrictions on the topology of Eq.(|^ domain used for Eq.(21) was violated. For example, 
when the shape of environment was rectangular and place held centers were located on rectangular lattice, 
grid held with a high degree of hexagonal symmetry were formed if sizes of place helds were not too big (less 
than 20 % of the size of environment) and coefficients near the second order terms of Eq.(|^ were sufficiently 
large (they were selected close to the transition to unimodal pattern formation mode) (supplementary Fig.SI) 

Restrictions on the maximal value of second order terms can be weakened by introducing additional 
plasticity rules to restrict maximal or minimal weights of individual synapses. For example, if the following 
parameter values are set: cr = 1, £2 = 2500, ri 2 + = 0, £1 = 0, 7y_|_ = 1, ? 7 _ = 1 , Ur = 0.13L in Eq. (21) and 
synaptic weights are additionaly limited from above by setting Wi = 0 for Wi > Wmax, Wmax = 0.005, then 
maximal ratio of second order term of Eq. (21) to the hrst order terms at the end of learning is approximately 
equal to 1. The resulting grid helds are shown in Fig. Hb- It is interesting to note that introduction of limits 
on maximal or minimal weights into Eq.(21) lacking second order terms (£2 = 0, 772 - 1 - = 0) is sufficient 
to obtain the model which could learn grid-like helds. Examples of grid helds obtained for such models 
are shown in Fig. [^-D. Parameter sets used in these simulations were the same as for Fig. but with 
£2 = 0 and Wmin = —0.001, cr = 0 for Fig. [^, £2 = 0 and Wmax = 0.005 for Fig. and £2 = 0, cr = 0 , 
Wmin = -0.001, Wmax = 0.005 for Fig. [^. 

On the contrary to the low sensitivity of grid held formation process to the second order terms coefficients, 
the ratio of hrst order terms coefficients could vary approximately within the limits of 77 _ = 0 . 6 ? 7 + and 
r]_ = 2 ? 7 _|_. But these limits correspond to physiologically different situations - Hebbian plasticity with 
large FTP component and very small LTD component for small rj- values and pure LTD with nonlinear 
dependence on presynaptic activity for large 77 - values. Examples of plasticity rate curves and grid helds 
obtained for different ratios of 77 - to 77 + are shown in Fig.[^ Here in Fig. dependencies of plasticity rate 
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Figure 2: Formation of grid fields in the models obtained from the basic model of neuron with associative plasticity by different 
modifications of terms describing local homeostatic plasticity (Eq. §)• Top row - Grid fields formed after 10000 time steps of 
simulation. Middle row - corresponding synaptic weights as a function of spatial location of presynaptic PCs fields. Lower row 
- Fouirer transforms of synaptic weights spatial distributions shown in the middle row. A. Imposing upper and lower bounds 
on the synaptic weights enable grid fields formation in the models with large second order terms. Parameters of the model 
(Eq .([^) were: a = 1, £2 = 2500, 772 + = 0, £1 = 0, 77 + = 1, rj— = 1 ^ ar = 0.13L. In addition, synaptic weights were limited 
from above by setting Wi = 0 for Wi > w-max, Wmax = 0.005. B. Imposing only lower bound on synaptic weights enable grid 
fields formation in models without other homeostatic plasticity terms. Parameters of the model (Eq. §) were the same as in A 
except for: cr = 0, £2 = 0, 772 + = 0, £1 = 0, 77 + = 1, 77_ = 1 , err = 0.13L, Wmin = —0.001. C. Models with upper limit imposed 
on synaptic weights learn inverted grid fields. Model parameters were the same as in A except for £2=0 and Wmax = 0.005. 
D. In the models with both upper and lower bounds imposed on synaptic weigths, clusters of synapses are formed in which 
size of grid fields and size of clusters can be regulated independently. Parameters of the model were the same as in A except 
for: £2 = 0, (T = 0 , Wjnin = —0.001, Wmax = 0.005 

of a PC-GC synapse on the distance of animal’s location from the center of corresponding place fields at 
fixed postsynaptic firing rate are shown for r]- = 0 . 6 ry+ and = 277 + by blue and green lines, respectively. 
It can be seen that in the case of small rj- the maximal value of negative plasticity rate (which could be 
considered as a result of LTD) is less than 7% of maximal positive plasticity rate. In the case of large ? 7 _ all 
plasticity changes could be considered as a result of LTD. Corresponding dependencies of wave amplitude 
change rates on the length of its wavevectors are shown in Fig.[^. Grid fields (top row), synaptic weigths as 
a function of spatial position of the centers of place fields (middle row) and wave amplitudes as a function of 
wavectors (lower row) are shown for the case rj- = 0 . 6 / 7 + Fig- [T foi' the case 77 - = 2rj+ in Fig. |T). 

It is possible to modify the model in such a way that formation of grid fields becomes possible due 
to interaction of LTD and BCM-like LTP terms without additional homeostatic plasticity restrictions on 
individual synaptic weights. This can be achieved, for example, if dependence of LTD rate on presynaptic 
activity is sublinear. Example of learning in such model is shown in Fig. . Learning equation was: 

Wi = -77_Ro7W + 772 +w^Q^w - cr ^ wlwi + £1 (42) 

k 

with parameters: a = 100, 772 + = 5000, ei = 5, rj- = 100 , Cr = 0.13L . First term of this equation was 
obtained by integrating over time the expression describing nonlinear LTD: 

« -77_RoiW (43) 

t 

where 

Ran = J ® L = Fo exp ^ IIri - f ^ (44) 

It should be noted that learning of grid fields in this model was possible only for narrow range of parameters. 
For example, the range of parameter Si values, for which the formation of grid fields was observed was 
between 3 and 5. 
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Figure 3: Formation of grid fields in the models with different ratios of potentiation and depression components of associative 
synaptic plasticity. Ratio of depression and potentiation components was regulated by changing parameter rj— from r]— = 0 . 677 + 
(blue lines in Fig.[3j\-,B and Fig. |3p ) to 77 - = 277 + (green lines in Fig.[3]4,B and Fig. [sJD) A. Dependencies of PC-GC synapse 
plastcity rate on the animal’s distance from the center of corresponding place fields at a fixed postsynaptic firing rate. It 
can be seen that in the case of small 77 — the maximal value of negative plasticity rate (which could be considered as those 
resulting from LTD) is less than 7% of maximal positive plasticity rate. In the case of large 77 — all plastic changes could 
be considered as resulting from LTD. B. Corresponding dependencies of wave amplitude change rates on the length of its 
wavevectors. C,D.Grid fieds (top row), synaptic weights as a functions of spatial position of the centers of place fields (middle 
row) and wave amplitudes as a function of wavectors (lower row). 


To confirm results obtained for simplified model described by Eq.(|^ in more realistic conditions we 
have modeled grid fields formation as a result of offline learning after environment exploration. In these 
simulations we assume that after exploration of diamond-shaped environment by an animal, episodic memory 
corresponding to the visited places could be spontaneously reactivated in such a way that relative activities 
of place cells during reactivation are similar to those at time of navigation. The physiologiacal basis for such 
reactivation could be for example a sharp wave activity observed in hippocampus during sleep and quiet 


wakefulness periods Csicsvari & Dupret (2014). 


We have observed that for sufflciently long path of movement, the process of spatially modulated fields 
formation was qualitatively similar to those described above for simplified model. For these simulations 
we directly solve equation obtained by simplification of Eqs. (§-§. After setting C(wt) = 

P{w^t,Xit) = BH{xit - 9)w% in Eqs. §-@ we obtain: 


Wit = -r]-Xztyt + m+Xitu'i + V+XuVt - O' X! '^ktWit + Bw%H {x^t - O) 


(45) 


where H is the Heaviside function. Offline learning was modeled by calculating synaptic weight changes at 
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Figure 4: Learning of grid fields in the model with sublinear LTD and BCM-like LTP terms without additional homeostatic 
plasticity restrictions on individual synaptic weigths. Panels are the same as in Fig. 


each iteration of the algorithm independetly for 5000 randomly selected places on the path of movement, 
and then updating synaptic weights using the sum of these weight changes. The set of parameter values 
was a = \, B = 1000, 772 + = 0, 77 + = 1, 77 - = 3, 0 = 0. Place fields were modeled with Gaussian functions 
Fi{v) « F exp(— 2 ^ Hr — ri||^) with F = 20. Place held size was varied from O.IL to 0.3L, where L = 1 
is size of diamond-shaped environment. Centers of place helds were initially selected at vertexes of 31x31 
hexagonal lattice which was embedded in environment. Then all place held centers were perturbed by 
addition of a random vector selected from Gaussian distribution with a = Q.bAL (AL - lattice spacing). 

Path of movement was obtained by the integration of the movement equations: 


rt+At =rt+ vt At 
Vi+At = 0.99vt -I- O.OlJvt, Svt Af (0,1) 


(46) 


At the borders of environment direction of motion was selected randomly and the speed ||vt|| was set to 0.1 
of its value just before collision with a border. Time step used for the integration of movement equations was 
At = 0.1. Average speed of motion along the path obtained after the solution of Eq.(46) was (||v||) = 0.06 
and the coefficient of variation of speed was CV = 0.64. 

Typical result of simulations is shown in Figj^ The total length of traveling path in a diamond-shaped 
environment of size L=1 was S=172. It can be seen that similarly to the case of model described by Eq.Q, 
formation of grid helds occurred in two stages: initially, a ring of waves with critical wavelength grew, 
then selection of six waves whose wavectors composed a hexagon took place. The difference with respect 
to simplihed model is that not all waves out of the hexagon of waves with maximal amplitudes dissapear 
completely at the end of learning and that the amplitudes of waves composing the hexagon are not equal. 
As a result, the gridness score of observed helds vary in a broad range and highly depends on grid size and 
path length. 

Learning of smaller size grid helds generally requires signihcantly longer paths. We observed that this 
requirement can be signihcantly weakened by introducing restrictions on a minimal weight of each synapse. 
We observed that in this modihed version of the model fast self-organization of grid helds could be achieved 
during online learning. 

To demonstrate this we have conducted a series of computational experiments. In these experiments 
the shape of environment was a square with a side L=I, place helds had nonperiodic Gaussian shape, and 
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Figure 5: Simulation of grid fields offline learning in the basic model of EC neuron resulting from memory reactivation after 
animal’s exploration of environment. Random path of the length S = 172L in a diamond-shaped labirinth of size L = 1 was 
simulated. Then patterns of PCs activity observed during navigation were reactivated in random order during the period of 
offline learning. Panels are the same as in Fig. [T] 


positions of their centers were selected randomly within the environment. Online learning was achieved by 
solving Eqs. 0-@ with global homeostatic plasticity switched off and local homeostatic plasticity in the 
form of lower bound on synaptic weights (i.e. equation P{wit,Xit) = —BH{wmin — Wit) was used instead of 
Eqs.Q). Online learning equation was: 

mt = Kt {-'q-Xityt + q+xftyt - BH {wmin - Wit)) (47) 

Parameter values were: 77 + = 1, ?]_ = 1.125, Wmin = —0.2, B = 1000. Decreasing learning rate was 
important for the fast formation of grid fields during online learning and at the same time for achievment of 
a high gridness score at the end of learning. Parameter Kt describes learning rate changes with time. It was 
selected in such a way that root mean square change of synaptic weights from a time moment t to t + At 
was equal to ^/{Awf} = \'^'q + 1.6 • 10“"^. Time step used for the integration of online learning equation 

was At = 0.1. Place fields were modeled as: Fi{r) « F exp(—Ik ~ with F = 1. 

Example of grid fields formation during online learning experiment is shown in Fig.[^ 

We found that in these experiments probability to obtain fields of size > 0.2L with gridness score higher 
than 0.5 at the end of online learning was high (in 69 cases of 100 experiments, gridness score in these 69 
experiments was Mean ± SD — 1.19 ± 0.3). Most part of gridness score growth (67%) in simulations with 
high score (> 0.5) at the end of learning was achieved during first 10 minutes of motion (with average speed 
V = 0.06L/s). The dependence of gridness score on the traveling path is shown in Fig. [Zb- 



Figure 6: Simulation of grid fields online learning in the model of EC neuron with bounded synaptic weights. Random path of 
the length S = 242L in a square labirinth of size L = 1 was simulated. The GC with 1000 synapses from PCs with place fields 
centers randomly distribured in environment was modeled. Panels are the same as in Fig. 
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Figure 7: Statistics of the quality of grid fieds observed in online learning simulations. A. Statistical box plots represent 
statistics of gridness scores for the grid fiel ds obtained in a series of N=100 online learning simulations for a model of EC 
neuron with bounded synaptic weights (Eq.( |47[ |). Each box corresponds to the animal’s movement paths of a different lengths 
embedded in the square environment of the size L = 1. On each box, the central red line coresponds to the median gridness 
score, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points not 
considered outliers, and outlier is indicated by the red cross. B. Examples of paths that were used in experiments described 
in A (represented by gray line, green circle - current position of the animal), corresponding to the mean path lengths 15 (top), 
30 (middle) and 80 (down). C Histograms represents distributions of gridness scores obtained in online learning experiments 
(N=100) corresponding to the mean path lengths 15 (left), 30 (middle) and 80 (right). 


Here, statistical boxes correspond to the results of simulations of online learning for N=100 different paths 
for 5 different path lengths. Examples of movement trajectories for path lengthes S=15 (upper panel),30 
(middle panel) and 80 (lower panel) are shown in Fig.[^. Distributions of gridness scores obtained in online 
learning experiments corresponding to the mean path lengths 15 (left), 30 (middle) and 80 (right) are shown 
in Fig. [7p. It can be seen from these experiments that large-sized fields > 0.2L, were Or- size of grid 
field) with relatively high positive gridness score (> 0.5) could be formed with high probability (19 of 100) 
after relatively short experience of navigation in a new environment (mean path length - S = 15L, mean 
time of motion with average speed v = 0.06L - T = 240 ) (Fig. [^, upper trace). Additional experiments 
demonstrated that formation of grid fields with relatively high gridness score was observed approximately 
for the same range of LTD/LTP ratios as those for which grid fields formation was possible in the basic 
model (Eq[^in the twisted torus domain. Sizes of grid fields formed varied from 11 for LTD/LTP=2.8 to 
23 for LTD/LTP=0.8 for the rectangular environment of size L=50 (Fig.S2 and Fig.S3). 

4. Discussion 

In this work we have described formation of grid-like spatially selective sensory fields in the group of 
models of EC neuron having linear activation function and recieving spatially modulated input from place 
cells or from neurons of different sensory systems. 
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We have shown that learning of hexagonaly symmetric fields in these models is a result of complex action 
of several simple biologicaly motivated synaptic plasticity rules. These rules include rules for associative 
synaptic plasticity similar to those widely used for description of linear LTD and LTP rule with supralinear 
dependence on presynaptic activity, or a combination of BCM-like rule for LTP with sublinear LTD rule. 
We observed that even solely LTD rule with nonlinear dependence of LTD rate on presynaptic activity, 
or a rule with very small LTD in comparison with LTP component, is sufficient for robust grid-like fields 
formation. 

From the perspective of pattern formation theory, the role in learning of the rules for associative synaptic 
plasticity with linear dependence on synaptic weigths is selection of patterns composed of a group of waves 
with similar wavelengthes. This critical wavelength determines the size of grid fields. 

Additional plasticity rules with nonlinear dependence on synaptic weights are required for selection of sub¬ 
group of six waves, whose wavevectors compose a hexagon, from the group of waves with critical wavelength. 
These rules can be a consequence of restrictions imposed on maximal or minimal synaptic weigths, of a small 
supralinearity of the dependence of synaptic plasticity rate on the synaptic weight, or nonlinearity of a neuron 
activation rule. BCM-like LTP term nonlinearly depends on synaptic weigths and as a result it can, together 
with sublinear LTD term, play a role in both selection of waves with critical wavelength and consequent 
formation of hexagonally symmetric pattern (Eq. 421. 

In addition to these results we have shown that within the basic model (Eq 45 and Eq[^ grid fields 
formation could be observed in experiments with a simulation of random movement in the environment, 
but learning in this model was slow because of restriction of the maximal values of second order terms of 
learning equations. 

On the other hand, the grid fields self-organisation model composed of LTD-like and LTP-like linear, with 
respect to synaptic weights, terms, and additional nonlinear local homeostatic plasticity term, which restrict 
minimal synaptic weight (Eq. @) , are capable of fast online learning of grid fields after relatively short 
expirience of navigation in a novel environment. ). The quality of grid fields was not lower than in offline 
learning experiments (for the same path length). The main difference between offline and online learning is 
that the latter requires selection of an appropriate schedule for learning rate changes at learning rates for 
which fast self-organization could be observed grid fields do not stabilize with a time. For online learning 
experiments described in the manuscript learning rate with a hyperbolic dependence on time was used (with 
ratio 10:1 between initial and final values). We found that main part of online learning could occur rapidly 
within 10-20 minutes of navigation in Im labyrinth with average speed 6cm/s (for 20cm fields ). Learning 
rate was the only parameter which should be additionally tuned for online learning. Nevertheless we assume 
that offline learning scenarios could be considered as a model of interaction of hippocampus with entorhinal 
cortex during hippocampal sharp-wave activity Csicsvari & Dupret (20141. 

The model described in this work differs from other published models of grid fields formation that are 
based on synaptic plasticity in PC-GC synapses. First, it does not require additional firing rate adaptation 


process as it is supposed to be important for grid fields learning in the model of Korpff and Treves (Kropff & 


Treves 2008 Si et al., 20121. As a result, learning in the proposed model is not sensitive to the inhomogeneity 


and anisotropy of the animal’s motion velocity distribution. In addition, the model allows learning on the 
basis of random reactivation of episodic memory of navigational information without requirment for recall 
of continuos fragments of movement trajectory. 

The speed of the grid field formation in modications of the proposed model, in which the lower boundary 


on synaptic weights is added (Eq.(47l), is similar to those observed for the model described in the work of 
Castro and Aguiar (Castro & Aguiar [2014), but learning rules used in our model have more general form 


and are better biologically motivated. 

Recently proposed model of Widloski and Fiete Widloski & Fiete (2014) significantly differs from the 
model proposed in this work. First, it assumes that EC neurons have unimodal place fields pre-existing 
before learning begins, and only the lateral EC synapses are modified during learning. Second, STDP rules 
are used in their model and different rules are required for different subpopulations of synapses. Mexican hat 
type lateral synaptic connections are formed as a result of learning and grid-like activity in their model is a 
result of a multi-bump attractor existing in the model after learning. It is interesting to note that frequency- 
dependent synaptic plasticity rule used in our work could be used instead of STDP rule within Widloski 
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and Fiete Widloski fc Fiete] ( |2014 ) model for learning of Mexican hat connectivity, and, if augmented with 
additional adaptation mechanism - for bump control learning. 

Some of the models (Eqj^ described in this work allowed detailed mathematical analysis. This property 
could be useful for the analysis of possible self-organization mechanism of modular structure of grid cell 
networks observed in experiments. These structures assume that several specific patterns of GCs network 
activity, including alignment of grid field orientations, should arise in a network of GCs. Using theoretical 
analysis we have demonstrated that models of GCs network composed of cells with plastic PC-GC connec¬ 
tions and fixed weak GC-GC connections could align orientations of their grid fields. This result is similar 


to those described in the works (Kropff & Treves, 2008 Si et al. 20121. 


From the experimental observations it is known that distribution of sizes of grid fields observed in MEG 
has discrete spectrum (Stensola et al. 2012). This property could have a simle explanation on the basis of 
the theory proposed in our work. It can be seen that the spatial distribution of synaptic weigths, which 
forms in the model with lower-bounded synaptic weights (Eq |47[ ) , has a spectrum consisting of peaks which 
form the hexagonal lattice (Eig.2B, middle panel). To reproduce discretness of grid field sizes, we propose a 
network consisting of a layer of neurons in EC (disconnected or weakly connected with each other) recieving 
strong synaptic input from a module of GCs with grid fields formed on the basis of model described by 
Eq.(47). It could be predicted that if synapses in this network have learning rule similar to PC-GC synapses 
rule, then GCs with a smaller grid fields sizes will be formed (with the ratio sqrt{3) : 1 or higher between 
sizes of fields). 

In conclusion we suggest that learning on the basis of simple and biologically plausible associative synap¬ 
tic plasticity rules could result in formation of grid-like fields in EC which receive synaptic inputs from 
hippocampal place cells and spatially modulated sensory inputs from different sensory systems. It is likely 
that such learning contributes to the formation of grid fields in early development and to the maintainence 
of normal grid cells activity patterns in familiar environments. 


Conclusions 

1. Learning of hexagonaly symmetric fields demonstrated in the model of EC neuron with linear activa¬ 
tion function. The neuron recieved spatially modulated input from place cells, or cells from sensory systems. 
Learning rules used in the model have a general form and are biologically motivated. 

2. Erom the perspective of pattern formation theory, associative synaptic plasticity, linear with respect to 
synaptic weigths, is required for the selection of patterns composed of waves with similar critcal wavelengthes. 
This critical wavelength determines the size of grid fields. Additional plasticity, which depends nonlinearly 
on synaptic weights, is required to make a selection of subgroup of six waves, whose wavevectors compose a 
hexagon, from the group of waves with critical wavelength. 

3. The model containing LTD-like and LTP-like linear, with respect to synaptic weights, terms, and 
additional nonlinear local homeostatic plasticity term, demonstarte fast online learning of grid fields in a 
novel environment even after a short navigational expirience. 

4. Aligment of grid field orientations in grid cell networks and discretness of grid fields sizes observed in 
experiments could be explained within the proposed model. 
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